"
RScript corresponding to the analyses in the article
  'The Gig Economy in Party Manifestos: Analysing the Salience of a New Issue 
  Across Europe'

This script is supplementary material to the above mentioned article. 
In combination with the data, it enables the full replication of all analyses.
If you encounter problems or find errors in the code, please contact the 
corresponding author of the article.
"


#### Analyses of Salience in Election Manifestos ####

##### Packages #####
if(!require(pacman)) install.packages("pacman")

library(pacman)
p_load(quanteda, quanteda.textmodels, quanteda.textstats, quanteda.textplots, 
       readtext, tidyverse, DT, stringr, ggplot2, sf, rnaturalearth, 
       rnaturalearthdata, lme4, readr, lmerTest,lmeresampler, MuMIn, sjPlot, car,
       corrplot, scales, psych, margins, HLMdiag)


#### Dictionary Analaysis ####
##### Data import and preprocessing #####
english_txts <- readtext::readtext("path_to_folder_with_txt-files/*",
                                   docvarsfrom = "filenames", 
                                   docvarnames = c("country", "party", 
                                                   "language"), 
                                   encoding = "UTF-8")

english_txts_corpus <- corpus(english_txts)

dfm_english_corpus <- english_txts_corpus %>%
  tokens(remove_punct = TRUE, remove_numbers = TRUE) %>%
  tokens_tolower() %>%
  tokens_remove(stopwords("english")) %>%
  tokens_wordstem() %>%
  dfm()


##### Setup of dictionary #####
myDict <- dictionary(list(gig = c("gig", "crowdwork*", "clickwork*", "flexw*"), 
                          employmentstatus = c("bogus", "self-employment"), 
                          contractstatus = c("contractor", "zero-hour", 
                                             "task-based","on-call","on-demand"),
                          platform = c("platform"),
                          workplatforms = c("Upwork", "Uber", "Deliveroo", 
                                            "Foodora")))

dictionary_gig <- dfm_lookup(dfm_english_corpus, myDict)
dictionary_gig_df <- convert(dictionary_gig, to = "data.frame")


##### Keyword-in-context anaylsis #####
corpus_toks <- tokens(english_txts_corpus, remove_punct = TRUE)

###### Direct address ######

kw_gigeconomy <- kwic(corpus_toks, 
                      pattern = phrase("gig e*"), 
                      window = 10)
kw_gigeconomy
df_gigecononomy <- data.frame(kw_gigeconomy)

kw_gigwork <- kwic(corpus_toks, 
                   pattern = phrase("gig w*"),
                   window = 10)
kw_gigwork
df_gigwork <- data.frame(kw_gigwork)

kw_platformwork <- kwic(corpus_toks, 
                        pattern = phrase("platform work*"),
                        window = 10)
kw_platformwork
df_platformwork <- data.frame(kw_platformwork)

kw_platformeconomy <- kwic(corpus_toks, 
                           pattern = phrase("platform eco*"), 
                           window = 10)
kw_platformeconomy
df_platformeconomy <- data.frame(kw_platformeconomy)

kw_crowdwork <- kwic(corpus_toks, 
                     pattern = phrase(c("crowdw*", "crowd-w*", "crowd w*")), 
                     window = 10)
kw_crowdwork
df_crowdwork <- data.frame(kw_crowdwork)

kw_clickwork <- kwic(corpus_toks, 
                     pattern = phrase(c("clickw*", "click-w*", "click w*")), 
                     window = 10)
kw_clickwork
df_clickwork <- data.frame(kw_clickwork)

kw_contractor <- kwic(corpus_toks, 
                      pattern = phrase("contractor"), 
                      window = 10)
kw_contractor
df_contractor <- data.frame(kw_contractor)

kw_zerohour <- kwic(corpus_toks, 
                    pattern = phrase(c("zero-hour", "zero hour",  "0-hour")), 
                    window = 10)
kw_zerohour
df_zerohour <- data.frame(kw_zerohour)

kw_ondemand <- kwic(corpus_toks, 
                    pattern = phrase(c("on-demand work*",  "on demand work*", 
                                       "work on-demand", "work on demand")), 
                    window = 10)
kw_ondemand
df_ondemand <- data.frame(kw_ondemand)

kw_oncall <- kwic(corpus_toks, 
                  pattern = phrase(c("on-call work*", "on call work*", 
                                     "work on-call", "work on call", 
                                     "on call contr*")), 
                  window = 10)
kw_oncall
df_oncall <- data.frame(kw_oncall)

kw_gigplatforms <- kwic(corpus_toks, 
                        pattern = phrase(c("Upwork", "Uber", "Deliveroo", 
                                           "Foodora")), 
                        window = 10)
kw_gigplatforms
df_gigplatforms <- data.frame(kw_gigplatforms)


df_allhits <- rbind(df_clickwork, df_contractor, df_crowdwork, df_gigecononomy, 
                    df_gigplatforms, df_gigwork, df_oncall, df_ondemand, 
                    df_platformeconomy, df_platformwork, df_zerohour)



####### Plot the number of hits per country #######
# Load Europe country data
europe <- ne_countries(scale = "medium", 
                       continent = "Europe", 
                       returnclass = "sf")

# Define the country scores
country_scores <- data.frame(name = c("Austria", "Belgium", "Czechia", "Denmark",
                                      "Estonia", "Finland", "Germany","Greece",
                                      "Ireland", "Italy", "Latvia", "Lithuania",
                                      "Netherlands", "Norway", "Poland", 
                                      "Portugal", "Romania", "Slovakia", 
                                      "Slovenia", "Spain", "Sweden", 
                                      "Switzerland", "United Kingdom"),
                             score = c(18, 41, 2, 1, 1, 14, 12, 1, 7, 2, 0, 0, 
                                       14, 9, 1, 5, 0, 2, 2, 3, 0, 9, 14))

# Merge with spatial data
europe <- europe %>%
  left_join(country_scores, by = "name")

# Define map boundaries for Europe
x_limits <- c(-25, 35)  
y_limits <- c(34, 72) 

# Plot the  map
map_hits <- ggplot(data = europe) + 
  geom_sf(aes(fill = score), color = "black", size = 0.2) +
  scale_fill_gradient(low = "gray80", high = "gray10", na.value = "white", 
                      name = "Number of hits") +
  coord_sf(xlim = x_limits, ylim = y_limits, expand = FALSE) +  
  theme_minimal() +
  theme(legend.position = "right",
        panel.grid = element_blank(),   
        axis.text = element_blank(),    
        axis.ticks = element_blank(),  
        axis.title = element_blank())

# Export the plot
ggsave(filename = "hit_map_europe.jpg", 
       plot = map_hits,
       width = 20,
       height = 20,
       units = "cm",
       dpi = 600)


###### Indirect address ######
kw_selfemployment <- kwic(corpus_toks, 
                          pattern = phrase(c("false self*", "bogus self*")), 
                          window = 10)
kw_selfemployment
df_selfemployment <- data.frame(kw_selfemployment)




#### Multilevel Modelling ####

##### Data import #####
df_regression <- read_delim("path_to_csv-file", 
                            delim = ";", na = "NA")

# construction of dummy variable for liberal parties
df_regression <- df_regression %>%
  mutate(liberal = ifelse(partyfamily == 4,1,0))

# rescaling all IVs (range [0;1])
df_regression_scaled01 <- df_regression %>% 
  mutate(partyfamily = scales::rescale(partyfamily),
         lrgen_score = scales::rescale(lrgen_score),
         lr_econ = scales::rescale(lr_econ),
         galtan = scales::rescale(galtan),
         seats_parliament = scales::rescale(seats_parliament),
         effective_parties = scales::rescale(effective_parties),
         GDP_pc_PPP = scales::rescale(GDP_pc_PPP),
         GDP_change_st = scales::rescale(GDP_change_st),
         GDP_change_st_percent = scales::rescale(GDP_change_st_percent),
         GDP_change_lt = scales::rescale(GDP_change_lt),
         GDP_change_lt_percent = scales::rescale(GDP_change_lt_percent),
         workingage_pop = scales::rescale(workingage_pop),
         internet_penetration_rate = scales::rescale(internet_penetration_rate),
         foreignborn_pop = scales::rescale(foreignborn_pop),
         EPL = scales::rescale(EPL),
         EPL_temp_contracts = scales::rescale(EPL_temp_contracts),
         unemp_expen = scales::rescale(unemp_expen),
         total_ws = scales::rescale(total_ws),
         LMP = scales::rescale(LMP),
         unemp_rate = scales::rescale(unemp_rate),
         invol_temp_contracts = scales::rescale(invol_temp_contracts),
         minimum_wage = scales::rescale(minimum_wage),
         union_density = scales::rescale(union_density),
         collective_bargaining = scales::rescale(collective_bargaining))



##### Descriptives #####
# correlation matrix with all standardised macro-level IVs
df_cor <- df_regression_scaled01 %>%
  select(effective_parties, GDP_pc_PPP, GDP_change_st, GDP_change_st_percent, 
         GDP_change_lt, GDP_change_lt_percent, proportional_es, workingage_pop, 
         internet_penetration_rate, foreignborn_pop, EPL, EPL_temp_contracts, 
         unemp_expen, total_ws, LMP, unemp_rate, invol_temp_contracts, 
         minimum_wage, union_density, collective_bargaining, citizenship_index)

cor_matrix <- cor(df_cor, use = "complete.obs")
cor_matrix <- round(cor_matrix, digits = 3)
colnames(cor_matrix) <- c("Number of effective parties", "GDP per capita (PPP)",
                          "GDP short-term change score", 
                          "GDP short-term change in percent", 
                          "GDP long-term change score",
                          "GDP long-term change in percent",
                          "Proportional electoral system",
                          "Share of Working-age population",
                          "Internet penetration rate",
                          "Share of foreign-born population",
                          "Employment protection legislation",
                          "EPL (Temporary contracts)",
                          "Unemployment expenditure",
                          "Total welfare state expenditure",
                          "Labour market policies expenditure",
                          "Unemployment rate",
                          "Share of involuntary temporary contracts", 
                          "Minimum wage", "Union density", 
                          "Collective bargaining coverage", 
                          "Ease of obtaining citizenship")

rownames(cor_matrix) <- c("Number of effective parties", "GDP per capita (PPP)",
                          "GDP short-term change score", 
                          "GDP short-term change in percent", 
                          "GDP long-term change score",
                          "GDP long-term change in percent",
                          "Proportional electoral system",
                          "Share of Working-age population",
                          "Internet penetration rate",
                          "Share of foreign-born population",
                          "Employment protection legislation",
                          "EPL (Temporary contracts)",
                          "Unemployment expenditure",
                          "Total welfare state expenditure",
                          "Labour market policies expenditure",
                          "Unemployment rate",
                          "Share of involuntary temporary contracts", 
                          "Minimum wage", "Union density", 
                          "Collective bargaining coverage", 
                          "Ease of obtaining citizenship")

jpeg(filename = "correlation_matrix.jpeg", 
      width = 20,
      height = 20,
      units = "cm",
      res = 600)
cor_plot <- corrplot(cor_matrix, is.corr = T, type = "upper", order = "AOE", 
                     tl.col = "black", tl.cex = 0.7)
dev.off()


##### Multilevel models with varying intercept #####
###### Null model ######
M0 <- lmer(direct_address ~ 1 + (1|country), data = df_regression_scaled01)
summary(M0)

# ICC null model
summary(M0)
0.03232/(0.03232+0.17241)


###### Party-level models ######
# model with all party-level variables without controls
M1 <- lmer(direct_address ~ partyfamily + liberal + lrgen_score + lr_econ 
           + galtan + seats_parliament + new_reelected + (1|country), 
           data = df_regression_scaled01)
summary(M1)
r.squaredGLMM(M1)
vif(M1)

# ICC model 1
0.0333/(0.0333+0.1509)


# party-level model without partyfamily, with covid dummy
M1.1 <- lmer(direct_address ~ liberal + lrgen_score + new_reelected 
             + seats_parliament + covid_dummy + (1|country), 
             data = df_regression_scaled01)
summary(M1.1)
r.squaredGLMM(M1.1)
vif(M1.1)

# ICC model 1.1
0.03543/(0.03543+0.15200)


# party-level models with different left-right variables
M1.2 <- lmer(direct_address ~ liberal + lr_econ + galtan + new_reelected 
             + seats_parliament + covid_dummy + (1|country), 
             data = df_regression_scaled01)
summary(M1.2)
r.squaredGLMM(M1.2)
vif(M1.2)

# ICC model 1.2
0.03077/(0.03077+0.15243)


# party-level models 1.1 and 1.2 with party and voting system controls
M2.1 <- lmer(direct_address ~ liberal + lrgen_score + new_reelected 
             + seats_parliament + effective_parties + proportional_es 
             + covid_dummy + (1|country), 
           data = df_regression_scaled01)
summary(M2.1)
r.squaredGLMM(M2.1)
vif(M2.1)


M2.2 <- lmer(direct_address ~ liberal + lr_econ + galtan + new_reelected 
             + seats_parliament + effective_parties + proportional_es 
             + covid_dummy + (1|country), 
           data = df_regression_scaled01)
summary(M2.2)
r.squaredGLMM(M2.2)
vif(M2.2)

# ICC model 2.2
0.02876/(0.02876+0.15131)


###### Country-level models ######
# economic variables: GDP p.c., GDP change scores, unemployment rate
M3 <- lmer(direct_address ~ GDP_pc_PPP + unemp_rate + covid_dummy + (1|country),
           data = df_regression_scaled01)
summary(M3)
r.squaredGLMM(M3)
vif(M3)

# model only using GDP
M3.1 <- lmer(direct_address ~ GDP_pc_PPP + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M3.1)
r.squaredGLMM(M3.1)
vif(M3.1)

# ICC model 3.1
0.02257/(0.02257+0.17264)


# models with short-term GDP change scores
M3.2 <- lmer(direct_address ~ GDP_change_st + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M3.2)
r.squaredGLMM(M3.2)
vif(M3.2)


M3.3 <- lmer(direct_address ~ GDP_change_st_percent + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M3.3)
r.squaredGLMM(M3.3)
vif(M3.3)


# models with long-term GDP change scores
M3.4 <- lmer(direct_address ~ GDP_change_lt + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M3.4)
r.squaredGLMM(M3.4)
vif(M3.4)

# ICC model 3.4
0.02742/(0.02742+0.17265)


M3.5 <- lmer(direct_address ~ GDP_change_lt_percent + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M3.5)
r.squaredGLMM(M3.5)
vif(M3.5)


# model with GDP and GDP short- and long-term change scores 
M3.6 <- lmer(direct_address ~ GDP_pc_PPP +  GDP_change_st + GDP_change_lt 
             + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M3.6)
r.squaredGLMM(M3.6)
vif(M3.6)


# model only using unemployment rate
M3.7 <- lmer(direct_address ~ unemp_rate + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M3.7)
r.squaredGLMM(M3.7)
vif(M3.7)

# ICC model 3.7
0.03299/(0.03299+0.17169)


# legal variables: EPL, EPL temporary contracts, citizenship index
M4 <- lmer(direct_address ~ EPL + EPL_temp_contracts + citizenship_index 
           + covid_dummy + (1|country),
           data = df_regression_scaled01)
summary(M4)
r.squaredGLMM(M4)
vif(M4)


# model only using EPL (temporary contacts) and citizenship index
M4.1 <- lmer(direct_address ~ EPL_temp_contracts + citizenship_index 
             + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M4.1)
r.squaredGLMM(M4.1)
vif(M4.1)

# ICC model 4.1
0.03663/(0.03663+0.17076)


# model only using EPL (temporary contacts)
M4.2 <- lmer(direct_address ~ EPL_temp_contracts + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M4.2)
r.squaredGLMM(M4.2)
vif(M4.2)


# welfare variables: unemployment expenditure, total welfare expenditure, LMP
M5 <- lmer(direct_address ~ unemp_expen + total_ws + LMP + covid_dummy 
           + (1|country),
           data = df_regression_scaled01)
summary(M5)
r.squaredGLMM(M5)
vif(M5)

# ICC model 5
0.02079/(0.02079+0.17309)


# model only using unemployment expenditure and total welfare spending
M5.1 <- lmer(direct_address ~ unemp_expen + total_ws + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M5.1)
r.squaredGLMM(M5.1)
vif(M5.1)


# model only using unemployment expenditure
M5.2 <- lmer(direct_address ~ unemp_expen + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M5.2)
r.squaredGLMM(M5.2)
vif(M5.2)


# collective action variables: union density, collective bargaining
M6 <- lmer(direct_address ~ union_density + collective_bargaining + covid_dummy
           +(1|country),
           data = df_regression_scaled01)
summary(M6)
r.squaredGLMM(M6)
vif(M6)


# model only trade union density
M6.1 <- lmer(direct_address ~ union_density + covid_dummy + (1|country),
             data = df_regression_scaled01)
summary(M6.1)
r.squaredGLMM(M6.1)
vif(M6.1)

# ICC model 6.1
0.04702/(0.04702+0.17039)


# country-level model with significant predictors and controls
M7 <- lmer(direct_address ~ GDP_pc_PPP + EPL_temp_contracts + unemp_expen 
           + covid_dummy + workingage_pop + internet_penetration_rate 
           + foreignborn_pop + (1|country),
           data = df_regression_scaled01)
summary(M7)
vif(M7)


M7.1 <- lmer(direct_address ~ GDP_change_lt + EPL_temp_contracts + unemp_expen 
           + covid_dummy + workingage_pop + internet_penetration_rate 
           + foreignborn_pop + (1|country),
           data = df_regression_scaled01)
summary(M7.1)
vif(M7.1)


M7.2 <- lmer(direct_address ~ EPL_temp_contracts + covid_dummy + workingage_pop 
             + foreignborn_pop + (1|country),
             data = df_regression_scaled01)
summary(M7.2)
vif(M7.2)



# models with significant party- and country-level variables
M8 <- lmer(direct_address ~ lrgen_score + seats_parliament + GDP_pc_PPP 
           + EPL_temp_contracts + unemp_expen + covid_dummy + (1|country),
           data = df_regression_scaled01)
summary(M8)
r.squaredGLMM(M8)
vif(M8)

# ICC model 8
0.02698/(0.02698+0.15234)


M8.1 <- lmer(direct_address ~ lr_econ + seats_parliament + GDP_pc_PPP 
           + EPL_temp_contracts + unemp_expen + covid_dummy + (1|country),
           data = df_regression_scaled01)
summary(M8.1)
r.squaredGLMM(M8.1)
vif(M8.1)

# ICC model 8.1
0.02488/(0.02488+0.14787)



##### Model diagnostics #####
l1_residuals_M1.2 <- hlm_resid(M1.2, level = 1, include.ls = F)
l1_residuals_M2.2 <- hlm_resid(M2.2, level = 1, include.ls = F)
l1_residuals_M3.1 <- hlm_resid(M3.1, level = 1)
l1_residuals_M4.1 <- hlm_resid(M4.1, level = 1)
l1_residuals_M5 <- hlm_resid(M5, level = 1)
l1_residuals_M6.1 <- hlm_resid(M6.1, level = 1)
l1_residuals_M8.1 <- hlm_resid(M8.1, standardize = T)
l1_residuals_M8.1_country <- hlm_resid(M8.1, level = "country", standardize = T)



##### Regression tabels #####
covariate.labels <- c('(Intercept)' = "Intercept",
                      liberal = "Dummy: Liberal party",
                      lr_econ = "CHES economic left-right placement",
                      galtan = "GAL/TAN placement",
                      new_reelected = "Dummy: Re- or newly elected party",
                      seats_parliament = "Number of seats in parliament",
                      effective_parties = "Number of effective parties",
                      proportional_es = "Dummy: Proportional electoral system",
                      GDP_pc_PPP = "GDP per capita (PPP)",
                      EPL_temp_contracts = "EPL (Temporary contracts)",
                      citizenship_index = "Ease of obtaining citizenship in country",
                      unemp_expen = "Unemployment expenditure",
                      total_ws = "Total welfare state expenditure",
                      LMP = "Labour market policies expenditure",
                      union_density = "Trade union density",
                      covid_dummy = "Dummy: Election during COVID-19 pandemic")


tab_model(M1.2, M2.2, M3.1, M4.1, M5, M6.1, M8.1,
          file = "regressiontable_1.html",
          p.style = "stars",
          emph.p = T,
          show.ci = F,
          digits = 3,
          show.se = T,
          linebreak = T,
          collapse.se = T,
          pred.labels = covariate.labels,
          dv.labels = c("M1", "M2", "M3", "M4", "M5", "M6", "M7"),
          order.terms = c(3,4,2,6,5,10,11,13,14,15,16,12,8,9,7,1))




##### Robustness checks #####
# logistic regressions
# M1.2
M1.2_log <- glmer(direct_address ~ liberal + lr_econ + galtan + new_reelected 
                  + seats_parliament + covid_dummy + (1|country), 
                  data = df_regression_scaled01,
                  family = binomial(link = "logit"))
ame_M1.2_log <- margins(M1.2_log)
summary(ame_M1.2_log)


# M2.2
M2.2_log <- glmer(direct_address ~ liberal + lr_econ + galtan + new_reelected 
                  + seats_parliament + effective_parties + proportional_es 
                  + covid_dummy + (1|country), 
                  data = df_regression_scaled01,
                  family = binomial(link = "logit"))
ame_M2.2_log <- margins(M2.2_log)
summary(ame_M2.2_log)


# M3.1
M3.1_log <- glmer(direct_address ~ GDP_pc_PPP + covid_dummy + (1|country),
                  data = df_regression_scaled01,
                  family = binomial(link = "logit"))
ame_M3.1_log <- margins(M3.1_log)
summary(ame_M3.1_log)


# M4.1
M4.1_log <- glmer(direct_address ~ EPL_temp_contracts + citizenship_index 
                  + covid_dummy + (1|country),
                  data = df_regression_scaled01,
                  family = binomial(link = "logit"))
ame_M4.1_log <- margins(M4.1_log)
summary(ame_M4.1_log)


# M5
M5_log <- glmer(direct_address ~ unemp_expen + total_ws + LMP + covid_dummy 
                + (1|country),
                data = df_regression_scaled01,
                family = binomial(link = "logit"))
ame_M5_log <- margins(M5_log)
summary(ame_M5_log)


# M6.1
M6.1_log <- glmer(direct_address ~ union_density + covid_dummy + (1|country),
                  data = df_regression_scaled01,
                  family = binomial(link = "logit"))
ame_M6.1_log <- margins(M6.1_log)
summary(ame_M6.1_log)


# M8.1
M8.1_log <- glmer(direct_address ~ lr_econ + seats_parliament + GDP_pc_PPP 
                + EPL_temp_contracts + unemp_expen + covid_dummy + (1|country),
                data = df_regression_scaled01,
                family = binomial(link = "logit"))
ame_M8.1_log <- margins(M8.1_log)
summary(ame_M8.1_log)



# regression table
tab_model(M1.2_log, M2.2_log, M3.1_log, M4.1_log, M5_log, M6.1_log, M8.1_log,
          file = "regressiontable_2.html",
          p.style = "stars",
          emph.p = T,
          show.ci = F,
          digits = 3, 
          show.se = T,
          linebreak = T,
          collapse.se = T,
          pred.labels = covariate.labels,
          dv.labels = c("M1", "M2", "M3", "M4", "M5", "M6", "M7"),
          order.terms = c(3,4,2,6,5,10,11,13,14,15,16,12,8,9,7,1))


# table of AMEs
ame_M1.2_df <- summary(ame_M1.2_log) %>% select(factor, AME=AME, SE=SE, p=p)
ame_M2.2_df <- summary(ame_M2.2_log) %>% select(factor, AME=AME, SE=SE, p=p)
ame_M3.1_df <- summary(ame_M3.1_log) %>% select(factor, AME=AME, SE=SE, p=p)
ame_M4.1_df <- summary(ame_M4.1_log) %>% select(factor, AME=AME, SE=SE, p=p)
ame_M5_df <- summary(ame_M5_log) %>% select(factor, AME=AME, SE=SE, p=p)
ame_M6.1_df <- summary(ame_M6.1_log) %>% select(factor, AME=AME, SE=SE, p=p)
ame_M8.1_df <- summary(ame_M8.1_log) %>% select(factor, AME=AME, SE=SE, p=p)




# bootstrapping
M1.2_boot <- bootstrap(M1.2, .f = fixef, type = "residual", B = 1000)
summary(M1.2_boot)
confint(M1.2_boot, type = "norm")


M2.2_boot <- bootstrap(M2.2, .f = fixef, type = "residual", B = 1000)
summary(M2.2_boot)
confint(M2.2_boot, type = "norm")


M3.1_boot <- bootstrap(M3.1, .f = fixef, type = "residual", B = 1000)
summary(M3.1_boot)
confint(M3.1_boot, type = "norm")


M4.1_boot <- bootstrap(M4.1, .f = fixef, type = "residual", B = 1000)
summary(M4.1_boot)
confint(M4.1_boot, type = "norm")


M5_boot <- bootstrap(M5, .f = fixef, type = "residual", B = 1000)
summary(M5_boot)
confint(M5_boot, type = "norm")


M6.1_boot <- bootstrap(M6.1, .f = fixef, type = "residual", B = 1000)
summary(M6.1_boot)
confint(M6.1_boot, type = "norm")


M8.1_boot <- bootstrap(M8.1, .f = fixef, type = "residual", B = 1000)
summary(M8.1_boot)
confint(M8.1_boot, type = "norm")



#### References: Packages in alphabetical order ####
citation("base")
citation("car")
citation("corrplot")
citation("DT")
citation("ggplot2")
citation("HLMdiag")
citation("lme4")
citation("lmeresampler")
citation("lmerTest")
citation("margins")
citation("MuMIn")
citation("rnaturalearth")
citation("rnaturalearthdata")
citation("psych")
citation("quanteda")
citation("quanteda.textmodels")
citation("quanteda.textplots")
citation("quanteda.textstats")
citation("readr")
citation("readtext")
citation("scales")
citation("sf")
citation("sjPlot")
citation("stringr")
citation("tidyverse")